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Mechanical response of an inclined frictional granular layer ap- 
proaching unjamming 
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Abstract - We present an orthotropic elastic analysis of frictional granular layers under gravity 
by studying their stress response to a localized overload at the layer surface for several substrate 
tilt angles. The distance to the unjamming transition is controlled by the tilt angle a with respect 
to the critical angle etc. We find that the shear modulus of the system decreases with a, but 
reaches a finite value as a — >■ Oc- We also analyze the vibration modes of the system and show 
that the soft modes play an increasing, though not crucial, role approaching the transition. 
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Various amorphous materials, and granular media in 
particular, exhibit a so-called jamming transition between 
rigid and flowing states. The nature of this transition 
has been investigated during the last decade, see recent 
reviews [l][2]. Most granular studies have focused on fric- 
tionless discs or spheres, typically controlled in volume 
fraction (j) or in pressure P [3H5j, showing that the jam- 
ming transition is critical (scaling exponents, diverging 
length scale) (3}[6J[7] and related to isostaticity (3|[8}{TT]. 
As the system looses its mechanical rigidity at the transi- 
tion, its shear modulus G is found to vanish as a power law 
with respect to the distance to jamming (jj — (j)c, where 0c 
is the critical volume fraction. Vibrational mode analysis 
have been performed on these frictionless granular systems 
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or similarly on soft glassy materials 
reporting an increasing number of soft modes, correspond- 
ing to collective low-energy motion of the particles, as the 
transition is approached. 

Much fewer studies have dealt with frictional grains 
in this context and they have mostly considered homo- 
geneous systems under isotropic pressure 11 18-22 . In 



the frictional case, the Liu-Nagel jamming concept 23 24 
must be revised 25 . In particular, jamming and isostatic 
points do not coincide any more [I], and one thus can ex- 
pect a finite shear modulus at the transition. In this letter. 



we consider static layers of frictional grains under gravity, 
by means of two-dimensional discrete element simulations 
(standard Molecular Dynamics [26"), and investigate their 
mechanical properties through the analysis of their stress 
response to a localized overload Fo at the layer surface. 
The layers are prepared at a fixed angle a with respect to 
the horizontal (see Fig. p^ for notations), and unjamming 
is approached as a is close to Uc, the critical value above 
which static layers cannot be equilibrated at that angle 
and always flow. Note that this unjamming point ttc is 
close in spirit to the situation of a jammed solid sheared 
up to its yield-stress 27 . It is also close, but different, to 



progressively tilted granular layers, which eventually loose 
their mechanical stability, see e.g., 28 29]. In our simula- 



tions, the volume fraction in the layer is fairly uniform all 
through the layer depth. Its value depends very weakly 
on the inclination angle, varying from (p ~ 0.823 at a = 
to (/) ~ 0.816 at ac (these are for the GG preparation, 
see below). These values are always larger ~ though not 
much - than the critical value, estimated in our system at 
(j)c ~ 0.815 [21 30J3l]. Then, here, the control parameter 
for the jamming/unjamming transition is the sole angle a. 
This situation is therefore qualitatively different to the ho- 
mogeneous configurations submitted to isotropic pressure 
cited above, and is effectively closer to an experimental 
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Fig. 1: (a) Normal stress response profile cr„„ for two values of a, and two values of 6 (see legend). Symbols: numerical data 
from the MD simulations (GG preparation, [32]). Solid lines: best elastic fit (see parameters in legend). The stresses are 
normalized by Fo/h and the distance by h. (b) Same for the shear stress response atn. 
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Fig. 2: Schematics of the system and notations, x is the hor- 
izontal axis, z is the vertical one, along which acts gravity g. 
The granular layer, of average thickness h, is inclined at an 
angle a with respect to horizontal, t and n are the axis respec- 
tively tangential and normal to the layer. A localized force 
Fq, which makes an angle 9 with respect to n, is applied on 
a grain close to the surface of the layer. The stress responses 
o"„„ and atn to this overload are measured at the bottom of the 
layer (fixed grains in blue). Axis (1,2), making an angle r with 
respect to {n,t), are those of the orthotropic elastic analysis. 



set-up. 

The numerical model is that described in 
N = 3600 polydisperse frictional discs coupled, when over- 



32 33 , with 



lapping, by normal and tangential linear springs, tangen- 
tial forces being limited by the Coulomb condition with a 
friction coefficient /i — 0.5. Two system preparations have 
been considered: a grain-by-grain (GG) and a rain-like 
(RL) deposition of the particles on a rough substrate con- 
sisting of fixed but size-distributed particles, inclined at 
the desired angle a. The layer is prepared when all grains 



of these protocols). No external pressure applied to the 
topmost layer of particles, i.e. the pressure in the system 
is due solely to the gravitational force acting on the parti- 
cles themselves. The typical thickness of the layer is ~ 23 
grain diameters, i.e. a system aspect ratio around 1/7. 
Above a certain inclination ac, the simulations preparing 
the packing before the response procedure do not converge 
towards a static layer - the grains always flow. The 'solid- 
liquid' transition is rather abrupt (Aa ~ 0.5°), allowing us 
for a well defined value of this critical angle: ac — 20.75° 
for the GG and ac — 20° for the RL preparations respec- 
tively. For the sake of concision, except otherwise stated, 
the displayed data are those obtained with the GG prepa- 
ration. 

Experimental and numerical works have shown that 
the linear stress response of granular systems to a point 
force is well described by (possibly anisotropic) elastic- 
ity 32 35-39 . Here we show that an orthotropic elastic 
model is able to reproduce quantitatively the normal and 
tangential stress bottom profiles, and that the correspond- 
ing shear modulus G decreases with a, in parallel to a de- 
crease of the coordination number and minor changes in 
the micro-structure (fabric). However, it does not vanish 
at the transition, but reaches a finite value as a —>■ etc. 
We finally analyze the vibration modes of the system and 
show that the soft modes play an increasing, though not 
crucial, role approaching the transition. 

Stress response functions. — Once a layer is de- 
posited, stabilized in an equilibrium state, an additional 
force Fq is applied on a grain close to the free surface. 



have reached static equilibrium (see 32 for a description 



and a new equilibrium state is reached (see 32 for de- 
tails). Taking the difference between the states after and 
before the overload, one can compute the contact forces 
in response to Fq. Introducing a coarse graining length 
w, the corresponding stress response can be determined. 
Taking w of the order of few mean grain diameters (here 
w = 6 (d)) as well as an ensemble averaging of the data 
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Fig. 3: Modulus ratio G/E\ determined from tlie elastic fit, 
as a function of a, for two different layer preparations (see 
legend). Inset: Same for E2/E\. 



(here N^ ~ 150 realisations for each tilt angle a), make the 
stress profiles quantitatively comparable to a continuum 
theory 
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such as elasticity, as discussed below. The 
amplitude of the overload was kept constant for all sim- 
ulations: Fq = 1.0 (m) 5, where (m) is the average mass 
of the grains. This value is sufficiently small to ensure a 
linear [40|[4T] and reversible response of the system for all 
values of a, including close to ac- 

Three examples of stress bottom profiles cr„„(i) and 
(Jtn{t) for the GG preparation are displayed in Fig. [l]for 
two values of the slope a of the layer and two values of 
the angle 9 that the overload force makes with the normal 
direction (see Fig. [2]). The stress profiles computed for the 
RL preparation have approximatively the same qualitative 
behavior. 

Orthotropic elastic analysis. — In order to repro- 
duce the stress response profile, isotropic elasticity is not 
enough (32], and we consider here orthotropic elasticity, 
characterized by a stiff axis (1) and a soft one (2) associ- 
ated to two Young moduli Ei and E2 < Ei. Two Poisson 
coefficients 1^12 and 1^21 rnust also be introduced. However, 
they are not independent and verify P12/E1 = i'2i/£'2- 
Finally, G is the shear modulus, so that the stress-strain 
relation, in the orthotropic axis writes: 



(1) 



Elastic energy is well defined if all moduli Ei,E2,G are 
positive and 1 — i^i2i'2i > 0. A last parameter of this 
modeling is the angle r that the axis (1,2) make with 
(n,i) (see Fig. [2]). 

The stress responses have been analytically computed 
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thickness, a numerical integration, in the spirit of |35| , 
must be done. Rough bottom boundary conditions (zero 
displacement) are imposed. We can adjust four dimension- 
less parameters (two modulus ratios G/Ei and E2/E1, a 
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Fig. 4: (a) Evolution of the coordination number Z (•) with 
the inclination of the layer a. Taking into account 'rattlers' 
(see text) gives a modified coordination number Z* (■). Right 
j/-axis: relative importance of friction mobilisation at contact 
(A), (b) Contact angle polar distributions (GG preparation) 
at a = 0,10,20°. Solid black line: fourth-order Fourier fit. 
Gravity is vertical (black arrow), (c) Fitted orthotropic elastic 
angle r as a function of a (*). The four characteristic angles 
of the contact angle distribution, computed with respect to the 
direction n, are also shown - these angles corresponds to the 
directions of the lobes, and those in between the lobes, see 
sketch and corresponding coloured arrows. 



Poisson coefficient 1^21 and the orthotropic angle r) to re- 
produce, for a given a, the numerical data for all 6. As 
shown in Fig. [I] (solid lines) , the fitting functions quanti- 
tatively describe the numerical data. 

In Fig. [3) we show the behavior of the elastic modulus 
ratios G/Ei and E2/E1 as the slope a is larger. They both 
decrease, almost by a factor of two, from a fairly constant 
value at small a to a lower one close to ac- In particular, 
G does not vanish close to the critical angle, in agreement 
with the observation that frictional granular systems re- 
main hyperstatic at the unjamming transition 
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Such a discontinuous behaviour has been seen in simula- 
tions by Otsuki and Hayakawa 
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43 for a semi-infinite layer (h — ^ cxd). For a finite layer in 25 



investigating the rhe- 
ology of sheared frictional grains close to jamming, and 
in experimentally created shear-jammed states reported 
The ratio G/E2 is almost constant (not shown). 



This behavior is also found for the RL preparation, al- 
though the variations of G/Ei are less pronounced. This 
softening is finally consistent with acoustic experiments on 
a granular packing in the vicinity of the transition l44| . 
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Microscopic variables. — In addition to the above 
global mechanical properties of the system, we have stud- 
ied the evolution of various microscopic quantities with a. 
The first one of interest is the coordination number Z , i.e. 
the average number of contacts per grain, here computed 
in the bulk of the layer, where it is fairly uniform - it ob- 
viously drops down close to the surface. Z monotonously 
decreases with a (Fig. [4k) and stays always far from the 
isostatic value Z{so = 3 (for frictional grains in 2D). As ev- 
idenced by the comparison of the curves in figures [3] and 
l4k, the modulus ratio G/Ei is not found to be a linear 
function oi Z — Ziso, in contrast with the finding of [18| 
on homogeneous frictional systems, close to isostaticity. 
Grains of the bulk that only carry their own weight do 
not contribute much to the global stability of the contact 
network. As for so-called rattlers in gravity-free packings 
(see 
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chap. 6), these grains can be removed from the 
contact counting, leading to a modified coordination num- 
ber of the layer Z* (see Fig. l4k). However, we have found 
that their number is roughly independent of a, so that 
the relation between G/Ei and Z* is not linear either. 
We have also studied the friction mobilisation at the con- 
tact level. In our MD simulations, the number of contacts 
with a ratio of the tangential force ft to the normal force 
/„ strictly equal to the microscopic friction ji is zero when 
static equilibrium is reached. However, some of them are 
effectively close to the Coulomb criterion, and we have 
computed the average (-^). This quantity, displayed in 
Fig. Hk, increases as a — > ac, but its overall variation is 
weak (see right y-scale). 

Finally, we have studied the contact angle distribution. 
Three of these distributions are represented as polar dia- 
grams for a = 0, 10 and 20 degrees in Fig. |4]d. The four 
strongly pronounced lobes are typical of the GG prepa- 
ration [34| (chap. 6). The vertical and horizontal direc- 
tions are always in between these lobes. When the layer 
is horizontal, the orthotropic stiff and soft directions are 
also found to be (almost) along the horizontal and ver- 
tical axis respectively. Note that the fitting procedure 
effectively gives here r — 93° in this case, while r = 90° 
(or 0°) would have been expected for symmetry reasons. 
This indicates the typical precision we have on the mea- 
sure of this orthotropic angle. Close to the critical slope, 
however, the orthotropic orientations are close to those of 
the lobes, the stiff one being in the direction of the slope. 
As evidenced in Fig|4j:, the transition between these two 
microscopic configurations occurs around a ~ 9°, i.e. well 
below ttc, in correspondence with the drop by a factor of 
2 of G/Ei between 8 and 10° (see Fig.|3]). 

Vibration mode analysis. — For each tilt angle 
a, we have computed the 3A'^ vibration modes of the 
layer, performing a harmonic approximation of the energy 
around the equilibrium point reached by deposition of the 
grains. We only studied layers obtained by GG deposi- 
tion. They do not show contacts being at the Coulomb 
threshold. This avoids the issue of the relevance of the 
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Fig. 5: Cumulative density of vibration modes Q{u)) for four 
different values of a (see legend). Inset: Number of modes njo, 
which are necessary to represent 70% of the displacement re- 
sponse, and corresponding highest frequency lDvo, as a function 
of Q. 



harmonic approximation for such contacts 19 . Each con- 



tact between two grains is thus being considered as two 
linear springs (normal and tangential stiffnesses fc„ and 
kt] here kt — 0.75fc„). The reference vibration frequency 
is that of a single spring-mass system wq = yjkn/ {m). 
We denote D{lu) = (37V)-i Emodc* '^(c:' - Co,) and Q{lu) = 
Jq D{ui') dcj', the density and cumulative density of states, 
where w = cj/wg is the normalized frequency. 

The function Q is displayed in Fig. [5] for four values of 
the slope a. No excess of vanishing frequencies are ob- 
served as the transition is approached, which is expected 
for a system that stays hyperstatic t42j. There is no sin- 
gle zero frequency mode appearing at the transition ei- 
ther. The criticality of the layer thus cannot be described 
by the emergence of an energetically costless direction in 
configuration space. However, some rather soft modes 
(0.1 < (D < 0.5) are enhanced when a — >■ ac- We have 
also computed the eigenfunctions corresponding to these 
modes. The displacement field in response to the overload 
with Fq can be decomposed on these functions, as they 
form a mathematical basis. We denote Ci the coefficients 
of this decomposition. A large part of the decomposition 
is owned by the 20 lowest frequency modes. Hence, we 
can compute the number of modes nyo which are necessary 
to represent 70% of the displacement field, and the high- 
est frequency wyp of these 7170 modes, by looking for the 
lowest frequency W70 such that X]i|iSi<is cf > 0.70. Both 
nyo and wyo are plotted against a in the inset of Fig. [5J 
While 7^70 stays approximatively constant and around 20, 
Wyo is found to decrease from ~ 0.3 by a factor of two 
when a varies from zero to ac, without however showing a 
clear singularity at the transition. This means that the re- 
sponse is progressively better represented by softer modes, 
even if these modes do not trigger the critical behavior of 
the layer close to the transition. 

Conclusions. — To sum up, we have simulated 2D 
frictional and polydisperse granular layers under gravity 
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inclined at an angle a, and investigated their mechanical 
and microscopic properties when the unjamming transi- 
tion is approached. This work tells us what to expect 
in real experiments, i.e. a layer that becomes elastically 
softer as a —>■ ac, but not to the point at which the system 
would loose its rigidity before avalanching. Progressive 
enhancement of rather soft vibration modes should be ob- 
servable, although they do not seem to play a particularly 
crucial role. 

An interesting continuation of this work is to study this 
system with smaller and larger values of the contact fric- 
tion coefficient /x. We have seen that GG layers prepared 
with fi = 0.1 have a significantly smaller critical angle 
ac — 14°, whereas those prepared with ^ = 10 {ac — 21°) 
do not differ much to the case fj, = 0.5. Similarly, the 
grains are more connected and more densely packed at 
that critical angle in layers with /i = 0.1 {Z ~ 3.6 and 
(f) ~ 0.833), whereas the values of Z and (p do not change 
much when the system is prepared with /i = 10. However, 
the systematic investigation of the role of fi on the me- 
chanical response of the granular layer, and on the elastic 
shear modulus in particular, is beyond the scope of the 
present paper. Another possible perspective could be to 
use granular simulations with a rolling resistance 45 in 
order to explore a wider range of (p, Z and a. 
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